// -*- mode: c++; indent-tabs-mode: nil; -*-
//
//The Biomolecule Toolkit (BTK) is a C++ library for use in the
//modeling, analysis, and design of biological macromolecules.
//Copyright (C) 2004-2006, Christopher Saunders <ctsa@users.sourceforge.net>,
//                         Tim Robertson <kid50@users.sourceforge.net>
//
//This program is free software; you can redistribute it and/or modify
//it under the terms of the GNU Lesser General Public License as published
//by the Free Software Foundation; either version 2.1 of the License, or (at
//your option) any later version.
//
//This program is distributed in the hope that it will be useful,  but
//WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//Lesser General Public License for more details.
//
//You should have received a copy of the GNU Lesser General Public License
//along with this program; if not, write to the Free Software Foundation,
//Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

#include <fstream>
#include <iostream>
#include <string>
#include <set>

#include <btk/core/algorithms/transforms.hpp>
#include <btk/core/atoms/pdb_atom.hpp>
#include <btk/core/molecules/polymer.hpp>
#include <btk/core/molecules/monomer.hpp>
#include <btk/core/io/pdb_system.hpp>

typedef BTK::ATOMS::PDBAtom atom;
typedef BTK::MOLECULES::Monomer<atom> monomer;
typedef BTK::MOLECULES::Polymer<monomer> polymer;
typedef BTK::IO::PDBSystem<polymer> pdb_system;

using BTK::ALGORITHMS::transform;

int main(int argc, char **argv) {

  if (argc < 3) {
    std::cerr << "usage: " << argv[0] 
         << " transform_file pdb_file [chain_id [chain_id]]\n";
    std::cerr << "\nMaps structure coordinates x to x' = Ux+T for\n"
         << "the rotation U and translation T specified in\n"
         << "transform_file.\n\n";
    return 1;
  }

  std::string transfile,pdbfile;
  std::set<char> chain_ids;

  transfile = argv[1];
  pdbfile = argv[2];

  if (argc >= 4) {
    for (int i = 3; i < argc; ++i) 
      chain_ids.insert(argv[i][0]);
  }

  pdb_system pdb(pdbfile);
  pdb_system::chain_iterator ci;

  if (!chain_ids.empty()) {
    ci = pdb.system_begin();

    while (ci != pdb.system_end()) {
      if (chain_ids.find(ci->chain_id()) == chain_ids.end()) {
        ci = pdb.erase(ci);
      } else ++ci;
    }
  }

  BTK::MATH::BTKMatrix U;
  BTK::MATH::BTKVector T;

  std::ifstream is(transfile.c_str());
  is >> U >> T;

  transform(pdb.structure_begin(),
            pdb.structure_end(),
            U,T);

  std::cout << pdb << std::endl;
}
